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Abstract 

Molecular dynamics simulations and integral equation calculations of a simple equimolar mixture 
of diatomic molecules and monomers interacting via attractive and repulsive short-range potentials 
show the existence of pattern formation (microheterogeneity), mostly due to depletion forces away 
from the demixing region. Effective site-site potentials extracted from the pair correlation func¬ 
tions using an inverse Monte Carlo approach and an integral equation inversion procedure exhibit 
the features characteristic of a short-range attractive and long-range repulsive potential. When 
charges are incorporated into the model, this becomes a coarse grained representation of a room 
temperature ionic liquid, and as expected, intermediate range order becomes more pronounced and 
stable. 
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I. INTRODUCTION 


Spontaneous pattern formation is a feature present in a diverse collection of physical, 
chemical and biological systems.- In spite of the diverse nature of these systems, the appear¬ 
ance of the emerging microphases is quite similar: in 2D systems circular droplets, stripes or 
“bubbles” occur, and in 3D systems one may hnd spherical droplets, sheets or tubes. In some 
cases the patterns appear as transient states due to energy or mass fluctuations that occur 
in the process of spinodal decomposition, but sometimes these states can be stabilized due 
to the presence of competitive interactions, in which one of the interactions is responsible 
for inhibiting the phase separation.-i^ 

The understanding of this self-organizing capability of soft and fluid matter is critical for 
a wide panoply of applications of great relevance nowadays. These self-assembly mechanisms 
play a crucial role in processes involving protein solutions in food products,-'^ therapeutic 
monoclonal antibodies,-^- nanolithography^ or gelation processes.— 

In the realm of colloidal science, systems with extremely short ranged repulsive inter¬ 
actions are often used as an experimental realization of the hard sphere fluid,— a system 
notorious for its theoretical interest. On the other hand, the addition of non-adsorbing poly¬ 
mers to the colloidal solution typically activates an attractive inter-particle interaction, due 
to the depletion mechanism. Moreover, changing the concentration and molecular weight 
of the polymer, the attraction range and strength of the colloid-colloid interaction can be 
tuned. Clustering is to be expected due to the presence of the attractive forces,— but 
in principle it would correspond to meta-stable states and/or irreversible processes of ki¬ 
netic nature. Nevertheless, microphases formed by clusters and percolating structures can 
be stabilized in protein solutions and colloid-polymer mixtures both in experiments^ and in 
theoretical descript ionsi^ due to the presence of additional repulsive interactions stemming 
from electrostatic forces. An extreme example of the stabilizing role of charges is the nanos- 
tructural organization that appears in room temperature ionic liquids (RTIL).— In fact, it 
has been shown, that long range repulsive interactions alone can give rise to nanostructural 
order,— the driving force of attractive interactions to induce spontaneous aggregation being 
replaced by external forces (e.g. pressure). 

In the case of colloidal systems, in which charged colloidal particles are screened by 
ions in the solvent, the colloid-colloid interaction has been shown on theoretical grounds to 
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be adequately represented by a Yukawa potentiaU^ii^ according to the Derjaguin-Landau- 
Verwey-Oberbeek (DLVO) theory. Following this, numerous works have resorted to poten¬ 
tials with a combination of short range attraction and a long range repulsion (SALR) in 
the form a double Yukawa,— or a Lennard-Jones (LJ) plus a Yukawa interaction^*^^*^^ in 
order to model the spontaneous emergence of microstructured patterns in fluids. On the 
other hand, back in 1999, Sear et al.,— made use of an empirical two exponential form with 
SALR characteristics in order to explain the experimental appearance of stable microphases 
of nanoparticles at the air-water interface. This potential has been studied in depth in model 
systems, both in bulk and in conhnement,— and as a rough approximation to account for 
vegetation patterns in ecosystems with limited resources.— 

In this work we will explore the possibility of pattern formation in a system in which 
only short ranged forces are present. Our model system, composed of heteronuclear dimers 
and monomers combines attractive and repulsive potentials, so as to mimic the interactions 
present in RTILs, but without electrostatic forces. To that aim we have performed extensive 
molecular dynamics simulations in the canonical and in the isothermal-isobaric ensembles. 
We will address the emergence of intermediate range order (IRO) analyzing the behavior 
of the partial, and concentration-concentration structure factors and performing a cluster 
analysis for various degrees of asymmetry in the sites of the diatomic particles. Reference 
Interaction Site Model (RISM) integral equation calculations have also been carried out, and 
are shown to agree remarkably well with the simulations results. By means of an Inverse 
Monte Carlo approach,— we have extracted effective interactions from the pair correlation 
functions of the simulated mixtures. For comparison, another set of effective potentials 
has been obtained from the RISM results using an integral equation inversion procedure. 
We will see, that despite the fact that all interactions at play are short ranged, their net 
effect leading to the pattern formation (microheterogeneity, or microstructure segregation 
at the nanoscale) translates into the appearance of effective interactions that agree with the 
characteristic trends of a short range attraction and a long range repulsion, i.e. a SALR 
potential. We have found that the effective potentials extracted from the simulation and 
those derived by the theoretical approach agree remarkably well. Finally, we have analyzed 
the role of charges on our model, which in fact by the addition of electrostatic site-site 
interactions becomes a rough representation of a RTIL. As expected, charges will be shown 
to enhance the pattern formation and the stability of the nanostructured phases. 
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The rest of the paper can be sketched as follows. In the next section we introdnce the 
model in full detail and briefly summarize the methodology. In Section III we introduce our 
most signihcant results. Conclusions and future prospects are to be found in Section IV. 

II. MODEL AND METHODS 

Our model consists in an equimolar fluid mixture of two different species, a two-site 
dimer AB, and a monomer C. The dimers are represented by a two center Lennard-Jones 
(LJ) site-site potential, in which the sites are separated by a distance 1. Our monomers also 
interact via LJ potentials. In all cases, the interactions are cut and shifted at a distance r^, 
by which the explicit form of the site-site potentials is 



( 1 ) 


and Uij{r) = 0 otherwise. Our model is to a certain degree inspired by the simple coarse¬ 
grained model for imidazolium based RTIL of Merlet et al..— We will see to what extent 


a simple model, with just two sites and purely short ranged interactions can reproduce 
the presence of nano-structural order as found in RTILs. To that aim we will however 
preserve the attractive/repulsive character of the interactions in the RTIL. In our model 
then, C monomers would correspond to anions, AB dimers to the molecular cations, with 
the imidazolium ring that contains the positive charge, being represented by site A, and 
the non-polar tail, by the larger site B. This implies that AA and CC interactions will be 
repulsive, BB and AC are attractive, finally BC and AB interactions are also repulsive. For 
the sizes of A and C particles we have chosen aAA = <^cc = 4 A, the elongation of the 
dimer I = sA. The AB distances of the dimers are hxed as constraints of the equations of 
motion. The LJ well is set to e = 2.092 kJ/mol, identical for all interactions. Since the 
size of the non-polar tail is essential to determine the nanostructural ordering,— we have 
considered various sizes for asB (with asB > cfaa always). For the attractive interactions 
we have truncated and shifted the LJ potential at Tc = SasB- For the repulsive interactions, 
we have simply used Tc = 2^/®crjj, thus defining purely repulsive soft spheres following the 
prescription of Weeks, Chandler and Andersen (WCA).— The complete set of parameters for 
all interactions is summarized in Table [B Finally, in order to analyze the effect of charges on 
the intermediate range order, we have considered explicitly the same model with a positive 


4 




charge +g on the A sites and a corresponding negative charge —q on the monomers. The 
value of q is varied between 0 and 0.25e, where e is the elementary electron charge. Again 
these values are of the same order as those considered in the model of Ref. 32 


A. Simulations and analysis 

We have carried out extensive molecular dynamics simulations of the system previously 
described using the LAMMPS package,—^— in the canonical and isothermal-isobaric ensem¬ 
bles using a Nose-Hoover thermostat and barostat.— Our samples contained 16384 parti¬ 
cles (samples of up to 65536 particles were investigated and no signihcant size dependence 
was found). For simplicity we considered equal masses for the three interaction centers: 
rriA = ms = rric = 16 g moR^. Initial thermalization runs at a temperature of 226 K were 
2 X 10® steps long, with a time step of 1 fs. Production runs were 5 x 10® steps long, and 
averages were carried out every 5000 steps. 

One of the problems one can encounter when performing canonical simulations in this type 
of system is the occurrence of phase transitions, either vapor-liquid equilibria or demixing. 
In order to guarantee that the states under consideration correspond to thermodynamic 
equilibrium conditions, and consequently any potential intermediate range order is not the 
result of a spinodal decomposition, we have run additional isothermal-isobaric simulations 
and analyzed the volume fluctuation of the samples. In this way one can avoid those states 
that lie inside the liquid-vapor spinodal. Moreover, one can compute the partial structure 
factors, dehned as 

Sij(k) = Xi6ij + XiXjp j {gij{r) - 1) e"*"''dr, (2) 

where p is the total number density, 6ij is a Kronecker 6, and Xi is the molar fraction of 
component i. Here sites A and B are considered as different particles and Pij is the atom- 
atom pair distribution function. Our samples are large enough to allow for an accurate 
integration of the pair distribution functions, and the results are consistent with direct k- 
sampling. Notice that as far as Eq.([2]) is concerned, xa = xb = xc = 1/3, hence in the 
large k limit all structure factors will tend to 1/3. From the partial structure factors it is 
possible to evaluate the concentration-concentration structure factor introduced by Bathia 
and Thornton,— for which we have dehned 

Scc{k) = x\j^Scc{k) + XcSAB-AB{k) - 2xABXcSc-AB{k), (3) 
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where now one has to consider explicitly the structure factors corresponding to the molecular 
species AB, and as a consequence xc = xab = 1/2. We can simply approximate Qab-ab = 
qbb and Qc-ab = Qcb, as if the scattering length or form factor of A sites was negligible 
compared to that of B sites. This is in principle not unreasonable given the much larger size 
of the B sites, but in a realistic situation one should take explicitly into account the true 
scattering lengths or form factors of sites A and B. Now one has to correct for the different 
values of the molar fraction when AB is considered as a single species and Eq. ([2]) is used 
in (j3]). In this way, hmfc_,.oo 5'cc(fc) = XcXab = 1/4. With all this in mind, the presence of a 
divergence when /c ^ 0 in Scc{k) is a signal of a demixing transition, so this quantity will 
be essential to assess the stability of the thermodynamic states chosen for our simulations. 

Finally, back to the vapor-liquid transition, one can analyze the corresponding k- 
dependent linear response susceptibility in density fluctuations, namely^ 


pksTxTik) 


|S(fc)| 

Ea6(^aXb)|S(fc)|af,’ 


( 4 ) 


whose k = 0 limit is precisely the isothermal compressibility. In Eq. (jl]) is Boltzmann’s 
constant, T the absolute temperature, and the elements of the matrix Sij are just the partial 
structure factors as dehned in Eq. ([2]). | ... | denotes the matrix determinant and | ... |a6 
the corresponding minor of the matrix S{k). The presence of a divergence -or a substantial 
increase in Xt( 0)- is a clear indication of the vicinity of a vapor-liquid transition. A careful 
monitoring of this quantity together with the use of NPT simulations provides a reliable 
assessment of the stability of the state points under consideration during the simulation 
runs. 

All systems and conditions studied in this work are summarized in Table [TTl In the case 
of system 8, when increasing the charge from O.lOe to 0.25e the conditions of temperature 
and density corresponding to systems 3, 6, and 7 lie in the two-phase region. Consequently 
we resorted to an isothermal-isobaric simulation at low positive pressure to achieve thermo¬ 
dynamic equilibrium conditions in our system with q = 0.25e. The final value of the total 
particle density achieved in this way is indicated in Table Ull 
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B. Inverse Monte Carlo method 


With the pair correlation functions produced along the simulation runs and the cor¬ 
responding statistical uncertainties calculated using block averages, we have used the In¬ 
verse Monte Carlo (IMG) procedure proposed by Almarza and Lomba^i in order to pro¬ 
duce single component site-site effective potentials able to reproduce the microscopic struc¬ 
ture exhibited by our mixture model. The procedure starts from a simple approximation 
(3ull^{r) = — log g{r) and proceeds to modify the pair potential along the simulation run in 
such a way that the calculated {r) matches the input g{r). Explicit details of the method 
can be found in Ref. 1^. In our case, we have used a total of 4000 particles. The procedure 
of inversion was carried out in 20 stages. In the last stages the effective potentials hardly 
varied, and the convergence between input and calculated ( 7 (r)’s according the prescription 


of Ref. 


31 


was achieved succesfully in all the cases. 


In this way, one can use as input of the IMC procedure either gAA^r), gBsi^), or gcc{'>^), 
and obtain a corresponding set of M^j^(r),and which will obviously be 

different, but in the case of emergence of intermediate range order should exhibit some 
common features. 


C. RISM integral equation 

The site-site correlations are obtained by solving the usual set of 2 equations, the site- 
site Ornstein-Zernike equation (SSOZ) and the closure equation, which we choose here to 
be the site-site hypernetted equation (SS-HNC). The SSOZ equation for the present system 
is explicitly given in the matrix form 


(W + |h)(W-‘ - |c) = I, (5) 

where the 3x3 matrix H (or C) has for elements Hij = hij{k){or Cij = Cij{k)), the Fourier 
transform (FT) of the site-site pair correlation functions hij{r) = gij{r) — 1 (or the direct 
correlation function Cjj(r)), where the index i,j stand for one of the sites A,B,C. The 
matrix W represents the intra-molecular correlations, which for the present system gives: 
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^ waa Wab Wac ^ ^ 


w = 


V 


/ 


1 jo{kl) 0 
jo{kl) 1 0 

0 0 1 




( 6 ) 


Wab wbb wbc 
Wac Wbc wcc 

where jo{x) is a spherical Bessel function. The matrix I is the identity matrix. The SS-HNC 
equations are written as 


gij{r) = exp 


Ujjjr) 

kBT 


+ hij{r) - Cij{r) , 


( 7 ) 


and there are 9 such independent equations to solve. 

Both equations are approximate ones, and their respective inconsistencies have been dis¬ 
cussed many times in the past literature.—"^ Based on empirical evidence from the literature, 
we expect that the correlations obtained through these equations for the present systems, 
both charged and uncharged, should be relatively good for the short range part, but perhaps 
not at long range. We are particularly interested to see if the correlations related to the 
appearance of the local structures can be reproduced by this theory. The structure factor 
dehned in Eq.(|3]) is the appropriate function for this purpose, as illustrated in the Results 
section. 

The practical solution of these equations consists in discretizing all the functions on an 
equidistant grid, both in r and k space. We use 2048 points with a r-grid of Ar = O.Ola^, 
which is enough for the present case to properly describe the asymptotic behavior of the 
correlations in direct and reciprocal space. The set of two equations are solved iteratively 
following techniques well documented in the literature. 

It is also possible to obtain the effective potentials which would correspond to the equiv¬ 
alent one-component representation of the system. This is achieved by imposing the pair 
correlation function to be the desired site-site correlation, namely g{r) = gxxif)-, in the set 
of the two integral equations for the 1-component system, and solve these equations for the 
direct correlation function and effective pair interaction. The direct correlation function can 
be obtained through the OZ equation for 1-component system (which is an exact relation); 


(1 -h psh{k)){l - psKk)) = 1, 


( 8 ) 


where h{r) = g{r) — 1 = hxxi.f') = gxx{i^) — 1, and the density ps is that of the effective 









1 component made solely of sites X. Once c(r) is obtained, one solves the HNC closure, 
which has the same form as Eq. Q, but now for the effective interaction Ueff{r) one gets: 

Ueff{r) = -ksT \kigxx{r) + hxx{r) - c(r)] (9) 


III. RESULTS 

A. Pair structure. 

Here we have analyzed the effect of the molecular geometry on the nanostructure for¬ 
mation changing the diameter of asB- We have first considered, asB = 8 A, 9 A, 10 A, 
and 12 A. Some snapshots of configurations for varying cbb are depicted in Figured! We 
have found that for asB > 9A clustering or microheterogeneity of C particles can only be 
appreciated when the packing of the B sites is so high that it resembles that of a solid. In 
fact in this case, the height of the hrst peak of SsBik) exceeds 2.7, which according to the 
Hansen-Verlet rule^ indicates that freezing conditions have been reached. Moreover, the 
prepeak in the structure factor characteristic of the presence of IRO is absent from SsBik). 
The clustering of C particles results from a merely steric effect, since these are restricted 
to occupy the holes between the large B particles. These effects can be appreciated in the 
snapshots of Figured! where the dense packing of B sites (red spheres) is readily apparent. 

For the reason mentioned above, we will concentrate on the results for asB = 8 A, and 9 
A . Already in the corresponding snapshot of Figured! one can appreciate the formation of a 
bicontinuous network of percolating clusters, connecting both AB dimers and C monomers. 
By bicontinuous network, we mean that the clusters formed by B-sites and C particles will 
be seen to both span practically the whole sample, forming two continuous interpenetrated 
percolating microphases. This can be analyzed from a more quantitative perspective by first 
taking a look at the corresponding pair distribution functions and partial structure factors, 
which are depicted in Figures d! and [3! respectively for Systems 1 to 6. Focusing first on the 
Qaa pair distribution function, one first appreciates the large exclusion hole after the hrst 
layer, which is a simple consequence of the large size of B-sites. Obviously the exclusion 
hole grows with the size of the B-sites, as can be seen when comparing Figures on the left 
and right columns. Correlations between A-sites extend up to hve (Taa^ and the width of the 
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gcc correlation is ~ 2acc- These features hint at the presence of some degree of IRO. B-B 
correlations (graphs in the middle row) behave like those of a dense fluid, and no apparent 
sign of clustering or IRO is evident. In contrast, the wide first peak of gcc is characteristic 
of clusters of particles confined in cavities, in this case formed by R-sites. This effect, as 
mentioned before is maximized for the largest asB- We will see later, that these clusters of 
partly occluded C-particles are connected, forming a three dimensional percolating structure. 

If we take now a look at the partial structure factors, we immediately appreciate a feature 
characteristic of the emergence of IRO, namely the presence of a prepeak at 0.25A“^. This 
corresponds to correlations in the range of 25A, the distance at which any sign of structure 
of the pair distribution function dies out. Interestingly, the prepeak is almost absent in 
Saa, except for a small maximum visible for the = 9 A and the highest density. This 
quantity shows otherwise very little structure for k > 0.5 A“^. As seen in the gAA^i the 
most relevant feature in the AA correlations is the exclusion hole due to the presence of the 
R-sites. In contrast, Sbb does exhibit a prepeak, even when no evidence of IRO was visible 
in gBB- This prepeak is more apparent in the monomer structure factor Sec- When the 
density is lowered the prepeak in the B-site structure factor shifts to lower Rvalues, and 
vanishes at p = O.OOlA”^. In the case of Sec-, the position of the prepeak is preserved, 
but its magnitude decreases. In Figure 0] the corresponding concentration-concentration 
structure factor is displayed. The prepeak at ko ~ 0.25A“^ is preserved, although its mag¬ 
nitude decreases when the total density is lowered. In contrast no increase when /c —)■ 0 is 
visible. This implies that we are encountering concentration fluctuations inducing spatial 
inhomogeneities, but no demixing transition. In Figure |5] we have plotted the k-dependent 
susceptibility corresponding to density fluctuations. The prepeak is visible except for the 
lowest density, which implies that density inhomogeneities with a spatial patterns are also 
correlated with the corresponding concentration inhomogeneities. But now, the /c —)■ 0 
behavior is different. As density is decreased the susceptibility (i.e. the isothermal com¬ 
pressibility) grows, an indication of the vicinity of a vapor-liquid transition. This means, 
that lowering the density from the value of p = O.OOlA”^ at the same temperature could 
move the system across the spinodal curve into the two-phase region. Our analysis indi¬ 
cates that the thermodynamic conditions we have simulated can be considered equilibrium 
states. Moreover, we have conhrmed that the results do not have a signiheant sample size 
dependence, by which metastability can also be ruled out. 
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The site-site correlation functions and structure factors obtained from the RISM theory 
are represented in dashed lines in Figs.2-3. It is seen that the agreement is excellent in most 
cases, particularly in what concerns the BB and CC correlations. The AA correlations are 
systematically underestimated near contact and overestimated at larger distances. The most 
significant differences are seen for the structure factors in Fig.3. Integral equations tend to 
exaggerate concentration fluctuations, and often tend to interpret small aggregate formations 
as such^i^ . We observe here a similar trend for the low density case p = O.OOlA”^, for which 
fluctuations compete the most with aggregate formation. The prediction of aggregation, 
through the pre-peak is in very good agreement with simulations for the highest density 
p = O.OOISA"^, precisely when the denser packing tends to favor aggregation. This is also 
in line with previous observations of similar type of behavior for model ionic liquids. These 
features are a direct consequence of the fact that the HNC closure approximation misses 
high order correlations, hence high order cluster contributions, which are represented in the 
bridge term 6p (r) that is neglected in the exponential of Eq. ([7]). We observe that in all 
cases, the k=0 behavior of the RISM structure factor always overestimates the concentration 
fluctuations. 


B. Effective pair potentials 

In Figure [6] we present the effective potentials obtained from the site-site pair distribution 
functions. By construction, using these effective potentials in a simulation for a single 
component system will lead to a pair distribution function coincident with the original site- 
site correlation of the mixture. This is one of the possible alternatives to reduce the behavior 
of a complex system to a simpler one component system. Other alternatives, such as the 
force-matching approach^ will lead to quantitatively different results, but certainly retaining 
the essential features of the effective potentials found here. Among these features, we see 
that in all cases the effective potential has a short range (extremely short in the case of AA 
potentials) attractive well and this is followed by a long range repulsive region, which extends 
to 20-30 A. The repulsive region of is much less visible and is illustrated in the inset. 
The repulsive range is more influenced by the change in the total density. The attractive part 
of AA and CC effective interactions is due to depletion forces (in this case the plain site-site 
interactions are repulsive). In the case of AA interactions, most of the attractive well is 
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masked by the excluded volume effect of the B sites in the AB molecules (the large repulsive 
potential between 5-15 A corresponds to the exclusion hole in qaa)- Note that even if in 
qbb long range correlations due to nanostructure organization are clearly not visible, there 
are long range repulsions in the BB effective potential, which are reflected in the prepeak in 
Sbb as an indication of IRO. The long range repulsion vanishes for p = 0.00lA“^, which we 
have seen is a state approaching the gas-liquid transition. 

Fig. [H] shows the effective pair potentials as obtained by the integral equation approach 
outlined in Section C. the comparison with the simulations is overall quite good in all cases. 
However, it is seen that the repulsive shoulder -which is the signature of the clustering 
ability- is always systematically underestimated by the theory. This is a direct consequence 
of the weaker tendency of the lET to predict clustering. 

Taken into account that B-sites are much larger that A-sites, we can think of our model 
as a system of B particles in a “sea” of C monomers, just like colloids in solution. Following 
Mani et ah— we can use a functional form of the type 


U{r)lkBT 


4ao 
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( 10 ) 


to represent the BB effective interactions. Note that given the large size of the B-sites, 
we have retained the repulsive part of the bare LJ interaction in order to account for the 
repulsive component of the effective potential. One can see that the fits of the effective 
interactions U^l/{kBT) to Eq. (ITOl) represented in Figure [7] are fairly accurate except for 
the minor inflection of the curve around 13 A . The parameters of the fit are collected in Table 
uni Notice that the exponent of the attractive LJ component, Oi deviates substantially from 
the standard value of 6 , being its range shorter as density increases. The range parameter 03 
grows considerably with the density, reflecting the increase of intermediate range ordering 
as the total density is increased. We observe that a single component representation of 
our system can be well performed by a standard SALR potential in which the long range 
repulsion has the form a Yukawa interaction, even when the original bare interactions in the 
mixture are relatively short ranged LJ potential. 
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C. Cluster analysis 


In order to go beyond the mere qualitative information provided by simulation snap¬ 
shots and the two-body level information furnished by pair distribution functions or site-site 
structure factors, we have also performed a geometric cluster analysis on the B sites and 
the C monomers, using different values for the link distance Vd- Essentially this distance 
dehnes two particles as linked, and in this work it has been dehned in terms of the position 
of the inflection point of the corresponding effective potentials depicted in Figure |6l We will 
use various values of in the range 10-12 A, for B-sites and C monomers, and 6-8 A for 
A-sites. The effects of the particular choice of on the cluster distribution will be analyzed. 
Specihcally, we have calculated the normalized cluster size distribution iV(s), as proposed by 
Stauffer.— This quantity is dehned as the fraction of particles contained in clusters of size s, 
i.e. N{s) = n(s)(s/7V), where n{s) is the number of clusters of size s. With this dehnition, 
^ A^(s) = 1. Of all the systems analyzed, in Figure [8] we have chosen to plot the results 
of System 6, which exhibits a signihcant prepeak in its partial structure factors. We ob¬ 
serve that the normalized cluster size distributions of both A, and, B-sites and C monomers 
present the same qualitative features: hrst one hnds a maximum for isolated particles which 
decays monotonously to zero at a value of cluster size, s, that strongly depends on Vd- This 
is a typical behavior of a non-associating huid, in which instantaneous clusters are created 
and destroyed as particles explore their conhgurational space. If stable hnite clusters were 
formed, the cluster size distribution should exhibit the corresponding maxima for the pre¬ 
ferred sizes. On the other end of the s-axis, interestingly one hnds large clusters that span 
all the simulation cell. Here N{s) shows little dependence on particularly for the B-sites 
and C monomers. Finally, the cluster size distribution of A and B sites is qualitatively very 
similar, which is understandable taking into account that both sites are linked into single 
molecular units. In the next section we will see that this symmetry is broken by the presence 
of charges and a new symmetry between A-sites and C particles emerges. 

Thus from our analysis, a more clear picture shows up, in which we have a large portion 
of the sample linked into microsegregated clusters forming bicontinuous structures, with a 
remnant of disconnected particles that form short lived structures up to tens or hundreds of 
particles depending on the choice of as one would expect in a non-associating fluid.. 
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D. The effect of charges 


Our previous results have shown that microheterogeneity, or stable intermediate range 
order can be induced by competing short range interactions in a simple mixture model of 
dimers and monomers. Our model was somehow inspired by a coarse grained representation 
of ionic liquids, which are in reality characterized by the presence of Coulombic interactions, 
absent from our model. An immediate question that deserves to be answered is then, 
how would the presence of charges affect the stability of the aforementioned bicontinuous 
structures ? To that aim we have carried out the corresponding analysis on systems 7 
and 8, that, as mentioned, correspond to system 3 with charges +q added to sites A and 
—q to the C monomers. For q = O.le, standard canonical molecular dynamics simulations 
were run. Recall that in the case of g = 0.25e, density had to be increased in order to 
move out of the vapor-liquid coexistence region. This was simply achieved by means of an 
isothermal-isobaric simulation run at the same T as the original system and a pressure of 
0.61 MPa, leading to a total p = 0.00195 A“^. In the snapshots of Figure [9] one can readily 
see that the charges enhance the formation of microstructural order, and particularly for the 
highest charge one see very well dehned stripes of C particles, stripes that now appear to 
be finite. A more clear picture emerges when taking a look at the partial structure factors, 
presented in left panels of Figure [101 Now the prepeak is perfectly dehned even for the 
Saa structure factor for the lowest charge, in contrast with the uncharged system Saa- The 
extremely large values of Sapiko) for 0.25 A“^, resemble Bragg peaks, and indicate 
the presence of quasi-periodic order in the microstructural domains. Moreover, if now one 
looks at the cluster size distributions plotted on the right panels of Figure HOl together 
with the percolating clusters, one hnds now a maximum centered at s ~ 20 for q = 0.25e 
for C and A-sites, which indicates the presence of hnite clusters of monomers and A-sites. 
This maximum is preserved in the results obtained for other charges up to g = 0.2e (not 
shown for the sake of brevity), to disappear for weaker Coulombic interactions. It is obvious 
that the net effect of charges on the microstructuring of our model mixture is to enhance 
the formation of nanostructures, also giving rise to the formation of finite size clusters for 
sufficiently high charges. In contrast, B-sites form a percolating bicontinuous structure 
coexisting with some disconnected B-sites or short lived aggregates. A-sites and C monomer 
form aggregates embedded in the percolating network of B-sites. All this suggests that the 
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network of B-sites forms cavities, with the A-sites pointing inside the cavity. This in turn 
is hlled by C monomers. This configuration is favored both by steric effects and by the net 
attraction between the positively charge A sites and negatively charged C monomers. 

On the other hand, despite the fact that A-sites form part of the AB dimers and C 
monomers are independent particles, due to the symmetry of the electrostatic interactions 
and the symmetry in shape and density -ctaa = crcc, Pa = Pc^, as the charges increase, AA 
and CC correlations become extremely similar -compare Saa and Sec in Figure fTOh. as one 
would encounter in a simple fully symmetric electrolyte. 

The next question is how this is all reflected on the effective potentials. These are plotted 
on FiguredH In all cases one observes the characteristic SALR structure, obviously being the 
CC and AA effective interactions those that are most affected by the introduction of charges. 
In spite of the fact that these two effective interactions result from the coarse graining of 
many body effects, the dominant role of electrostatic interactions already reflected in the 
partial structure factors leads to surprisingly similar effective potentials when charges are 
present. On the other hand the changes in are just quantitative. The attractive part is 
hardly influenced by the charges, since it results mostly from the depletion interactions and 
the bare attractive ubb- The long range repulsion is enhanced, and as the charge reaches 
q = 0.25e oscillations appear. These oscillations recall the Friedel oscillations characteristic 
of effective cation-cation potentials in liquid metals.— In the latter instance, the oscillations 
result from the quantum nature of the electrons. Here they result from the interplay of the 
Coulombic interactions and depletion forces. Thus for sufficiently large charges the long 
range attractive interaction between C, and A sites propagates through the AB bonds and 
induces the attraction well around 30A as a result of a many body effect. 


IV. CONCLUSIONS 

In summary, we have shown that a simple mixture of heteronuclear AB dimers and C 
monomers, with short range attractive and repulsive interactions designed so as to mimic 
the interactions present in RTILs, can give rise to the presence of nanostructural order 
in the form of micro-segregation in bicontinuous structures. This in turn translates into 
the characteristic presence of a prepeak in the site-site structure factors. These features 
are found both in simulation and in the integral equation results. The effective site-site 
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potentials extracted from the pair distribution functions by means of an IMC and integral 
equation approach, display the characteristic features of the SALR interactions, with the 
repulsive long range increasing as the total density (and hence the ordering) increases. The 
addition of charges to the model enhances the nanostructural order. When charges are large 
enough, one hnds well structured phases in which bicontinuous structures coexist with hnite 
size aggregates of monomers, caged in cavities formed by a network of the large uncharged 
sites, and with the cationic sites facing the inner part of the cavity. The effect of charges 
on our simple and rather symmetric model induces the symmetrization of the correlations 
of the anionic monomers and the cationic sites. The microscopic structure formed by the 
uncharged sites (apolar head in the RTILs) retains its bicontinuous nature and even if it is 
stabilized and enhanced by the charges is still mostly dominated by depletion effects and the 
bare short range attraction of the B-sites. In this regard, it is interesting to note that the 
appearance of a pre-peak in the wide angle scattering experiments and computer simulations 
of RTILS have been a subject of much investigations^*^ and has been related to the charge 
ordering and the subsequent appearance of segregated charged and uncharged molecular 
domains. Our work presents a unihed view of microsegregated bi-continuous domains, pre¬ 
peaks in structure factors and SARL type interactions, which are common to many complex 
systems. 

Obviously a much richer variety of structures would result from longer attractive un¬ 
charged tails, beyond the single B-site model used here. On the other hand, our simple 
model when reduced to two dimensions most likely will also give rise to more complex struc¬ 
tures, which in three dimensions are hindered by entropic effects. This is certainly a problem 
relevant to the behavior at interfaces which we intend to address in the future. 
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TABLE I. Lennard-Jones potential parameters. 


Particle i 

Particle j 

Interaction 

e (kJ/mol) 

(Tij 

Tc 

A 

A 

repulsive 

2.092 

4.0 A 


A 

B 

repulsive 

2.092 

{o-AA + 

2^/6 • uab 

A 

C 

attractive 

2.092 

4.0 A 

^ ■ O'BB 

B 

B 

attractive 

2.092 

O'BB 

^ ■ O'BB 

B 

C 

repulsive 

2.092 

{o'BB + 0-Cc)l‘^ 

2^/6 • OBC 

C 

C 

repulsive 

2.092 

4.0 A 

2^/6 • free 


TABLE 11. Potential parameters and thermodynamic state variables for the systems under study. 



Potential 

Thermodynamic state 



kl(e) 

(JB (A) p{k 3) 

T(K) 

P(MPa) 

System 1 

0 

8.0 

0.001 

226.4 

27.05 

System 2 

0 

8.0 

0.00125 

226.5 

39.5 

System 3 

0 

8.0 

0.0015 

226.5 

59.4 

System 4 

0 

9.0 

0.001 

226.5 

30.4 

System 5 

0 

9.0 

0.00125 

226.4 

53.2 

System 6 

0 

9.0 

0.0015 

226.4 

96.7 

System 7 

0.1 

8.0 

0.0015 

226.4 

39.4 

System 8 

0.25 

8.0 

0.00195 

226.3 

0.61 
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TABLE III. Parameters of tha SALR effective interaction (jlOl) between B sites fitted to the data 
extracted from the IMC procedure. Note that the potential is scaled with ksT, by which oq is 
dimensionless. 


ao ai 02 (A) 03 (A) 

System 1 1.788 8.185 3.749 4.297 
System 2 2.066 8.667 0.843 7.282 
System 3 3.578 9.927 0.231 14.816 



(a)fTBB = SA (b)crBB = I2A 

FIG. 1 . Snapshots of configurations for total particle density p = 0.00125A“^ and temperature 
T = 226.45AT for two B-site diameters. As the size of B-sites grows C monomers cluster in the 
cavities formed by the B-sites due to excluded volume effects. All other diameters and total density 
are kept fixed, uaa = <^cc = 4.0 A . 
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r /A r/A 

{a,)aBB = SA {h)cFBB = 9A 

FIG. 2. The figures show the radial distribution functions for A, B and C particles respectively. 
Column (a) corresponds to ctbb = SA for system 3 (theory vs. simulation) and column (b) presents 
the simulations results for systems 4 to 6 for aBB = 9A. Total density is indicated in the legend. 
Simulation results are represented by solid lines and dash-dotted curves correspond to integral 
equation calculations. 
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FIG. 3. The figures show the structure factors for A, B and C particles respectively. Column 
(a) corresponds to asB = SA for system 3 (theory vs. simulation) and column (b) presents the 
simulations results for systems 4 to 6 for = 9A. Total density is indicated in the legend. 
Simulation results are represented by solid lines and dash-dotted curves correspond to integral 
equation calculations. 
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k(1/A) 


FIG. 4. Concentration-concentration structure factor for the Systems 1, 2 and 3. 
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k (1/A) 


FIG. 5. Isothermal compressibility as a function of k for the Systems 1, 2 and 3. 
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r/A r /A 

= SA (b)(7BB = 9A 

FIG. 6. Effective potentials for A, B and C particles respectively. Column (a) corresponds to 
<7 _b_b = SA for system 3 (theory vs. simulation) and column (b) presents the simulations results 
for systems 4 to 6 for (Tbb = SA. Total density is indicated in the legend. Simulation results are 
represented by solid lines and dash-dotted curves correspond to integral equation calculations. 
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FIG. 7. B-B effective interaction for systems 1 to 3, fitted to a generalized LJ+Yukawa interaction 
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FIG. 8. Normalized cluster distribution function for A, and B sites, and C monomers of System 3. 
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{si)qA = O.lOe; qc = —O.lOe 



(b)gA = 0.25e; qc = -0.25e 


FIG. 9. Snapshots of the equimolar mixture of AB dimers and C monomers with embedded charges 


(indicated on the figures). 



■i 


(b) Cluster size distribution 



FIG. 10. (a) Charge dependence of the partial structure factors for A (top), B(middle) and C 
(bottom) particles (b) Charge dependence of the cluster size distribution. Charge magnitudes are 
specified in the legend. 
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FIG. 11. Charge dependency of the effective potentials for A (top), B(middle) and C (bottom) 
particles. Charge magnitudes are specified in the legend. Values of rd correspond to the inflexion 
points of the effective potentials in their first minimum, i.e., — A) = 6 A, rd{B — B) = llA, 

rd{C-C) = 10 A 
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